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ABSTRACT 

The early Solar System contained short-lived radionuclides such as ^'^Fe 
{ti/2 = 1-5 Myr) whose most likely source was a nearby supernova. Previous 
models of Solar System formation considered a supernova shock that triggered 
the collapse of the Sun's nascent molecular cloud. We advocate an alternative hy- 
pothesis, that the Solar System's protoplanctary disk had already formed when a 
very close (< 1 pc) supernova injected radioactive material directly into the disk. 
We conduct the first numerical simulations designed to answer two questions re- 
lated to this hypothesis: will the disk be destroyed by such a close supernova; 
and will any of the ejecta be mixed into the disk? Our simulations demonstrate 
that the disk does not absorb enough momentum from the shock to escape the 
protostar to which it is bound. Only low amounts (< 1%) of mass loss occur, due 
to stripping by Kelvin-Helmholtz instabilities across the top of the disk, which 
also mix into the disk about 1% of the intercepted ejecta. These low efficiencies 
of destruction and injectation are due to the fact that the high disk pressures 
prevent the ejecta from penetrating far into the disk before stalling. Injection 
of gas-phase ejecta is too inefficient to be consistent with the abundances of ra- 
dionuclides inferred from meteorites. On the other hand, the radionuclides found 
in meteorites would have condensed into dust grains in the supernova ejecta, and 
we argue that such grains will be injected directly into the disk with nearly 100% 
efficiency. The meteoritic abundances of the short-lived radionuclides such as 
^°Fe therefore are consistent with injection of grains condensed from the ejecta 
of a nearby (< 1 pc) supernova, into an already-formed protoplanctary disk. 

Subject headings: methods: numerical — shock waves — solar system: formation — stars: 
formation — supernovae: general 
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1. Introduction 

Many aspects of the formation of the Solar System are fundamentally affected by the 
Sun's stellar birth environment, but to this day the type of environment has not been well 
constrained. Did the Sun form in a quiescent molecular cloud like the Taurus molecular 
cloud in which many T Tauri stars are observed today? Or did the Sun form in the vicinity 
of massive O stars that ionized surrounding gas, creating an H ll region before exploding as 
core-collapse supernovae? Recent isotopic analyses of meteorites reveal that the early Solar 
System held live ^°Fe at moderately high abundances, ^°Fe/^^Fe ~ 3 — 7 x 10~^ (Tachibana 
& Huss 2003; Huss & Tachibana 2004; Mostefaoui et al. 2004, 2005; Quitte et al. 2005; 
Tachibana et al. 2006). Given these high initial abundances, the origin of this short-lived 
radionuclide (SLR), with a half-life of 1.5 Myr, is almost certainly a nearby supernova, and 
these meteoritic isotopic measurements severely constrain the Sun's birth environment. 

Since its discovery, the high initial abundance of ^°Fe in the early Solar System has 
been recognized as demanding an origin in a nearby stellar nucleosynthetic source, almost 
certainly a supernova (Jacobsen 2005; Goswami et al. 2005; Ouellette et al. 2005; Tachibana 
et al. 2006, Looney et al. 2006). Inheritance from the interstellar medium (ISM) can be 
ruled out: the average abundance of ^°Fe maintained by ongoing Galactic nucleosynthesis in 
supernovae and asymptotic-giant-branch (AGB) stars is estimated at ^^Fe/^^Fe — 3 x 10~^ 
(Wasserburg et al. 1998) to 3 x 10"''' (Harper 1996), lower than the meteoritic ratio. 
Moreover, this ^°Fe is injected into the hot phase of the ISM (Meyer & Clayton 2000), and 
incorporation into molecular clouds and solar systems takes ~ 10^ years or more (Meyer & 
Clayton 2000; Jacobsen 2005), by which time the ^°Fe has decayed. A late source is argued 
for (Jacobsen 2005; see also Harper 1996, Meyer & Clayton 2000). Production within the 
Solar System itself by irradiation of rocky material by solar energetic particles has been 
proposed for the origin of other SLRs (e.g., Lee et al. 1998; Gounelle et al. 2001), but 
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neutron-rich ^°Fe is produced in very low yields by this process. Predicted abundances are 
60pg^56pg ^ 10~^^, too low by orders of magnitude to explain the meteoritic abundance 
(Lee et al. 1998; Leya et al. 2003; Gounelle et al. 2006). The late source is therefore a stellar 
nucleosynthetic source, either a supernova or an AGB star. AGB stars are not associated 
with star-forming regions: Kastner & Myers (1994) used astronomical observations to 
estimate a firm upper limit of ^ 3 x 10^® per Myr to the probability that our Solar System 
was contaminated by material from an AGB star. The yields of ^°Fe from an AGB star also 
may not be sufficient to explain the meteoritic ratio (Tachibana et al. 2006). Supernovae, 
on the other hand, are commonly associated with star-forming regions, and a core-collapse 
supernova is by far the most plausible source of the Solar System's ^°Fe. 

Supernovae are naturally associated with star-forming regions because the typical 
lifetimes of the stars massive enough to explode as supernovae ( ^ 8 Mq) are ^ 10^ yr, too 
short a time for them to disperse away from the star-forming region they were born in. 
Low-mass (~ 1 Mq) stars are also born in such regions. In fact, astronomical observations 
indicate that the majority of low-mass stars form in association with massive stars. Lada & 
Lada (2003) conducted a census of protostars in deeply embedded clusters complete to 2 kpc 
and found that 70-90% of stars form in clusters with > 100 stars. Integration of the cluster 
initial mass function indicates that of all stars born in clusters of at least 100 members, 
about 70% will form in clusters with at least one star massive enough to supernova (Adams 
& Laughhn 2001; Hester & Desch 2005). Thus at least 50% of all low-mass stars form in 
association with a supernova, and it is reasonable to assume the Sun was one such star. 
Astronomical observations are consistent with, and the presence of ^"^Fe demands, formation 
of the Sun in association with at least one massive star that went supernova. 

While the case for a supernova is strong, constraining the proximity and the timing 
of the supernova is more difficult. The SLRs in meteorites provide some constraints on 
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the timing. The SLR ^°Fe must have made its way from the supernova to the Solar 
System in only a few half-lives; models in which multiple SLRs are injected by a single 
supernova provide a good match to meteoritic data only if the meteoritic components 
containing the SLRs formed < 1 Myr after the supernova (e.g., Meyer 2005, Looney et al. 
2007). The significance of this tight timing constraint is that the formation of the Solar 
System was somehow associated with the supernova. Cameron & Truran (1977) suggested 
that the formation of the Solar System was triggered by the shock wave from the same 
supernova that injected the SLRs, and subseqeuent numerical simulations show this is a 
viable mechanism, provided several parsecs of molecular gas lies between the supernova 
and the Solar System's cloud core, or else the supernova shock will shred the molecular 
cloud (Vanhala & Boss 2000, 2002). The likehhood of this initial condition has not yet been 
estabhshed by astronomical observations. Also in 1977, T. Gold proposed that the Solar 
System acquired its radionuclides from a nearby supernova, after its protoplanetary disk had 
already formed (Clayton 1977). Astronomical observations strongly support this scenario, 
especially since protoplanetary disks were directly imaged ~ 0.2 pc from the massive star 
9^ Ori C in the Orion Nebula (McCaughrean & O'Dell 1996). Further imaging has revealed 
protostars with disks near (< Ipc) massive stars in the Carina Nebula (Smith et al. 2003), 
NGC 6611 (Oliveira et al. 2005), and M17 and Pismis 24 (de Marco et al. 2006). This 
hypothesis, that the Solar System acquired SLRs from a supernova that occurred < 1 pc 
away, after the Sun's protoplanetary disk had formed, is the focus of this paper. 

In this paper we address two main questions pertinent to this model. First, are 
protoplanetary disks destroyed by the explosion of a supernova a fraction of a parsec away? 
Second, can supernova ejecta containing SLRs be mixed into the disk? These questions 
were analytically examined in a limited manner by Chevalier (2000). Here we present the 
first multidimensional numerical simulations of the interaction of supernova ejecta with 
protoplanetary disks. In §2 we describe the numerical code, Perseus, we have written to 
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study this problem. In §3 we discuss the results of one canonical case in particular, run 
at moderate spatial resolution. We examine closely the effects of our limited numerical 
resolution in §4, and show that we have achieved sufficient convergence to draw conclusions 
about the survivability of protoplanetary disks hit by supernova shocks. We conduct a 
parameter study, investigating the effects of supernova energy and distance and disk mass, 
as described in §5. Finally, wc summarize our results in §6, in which we conclude that 
disks are not destroyed by a nearby supernova, that gaseous ejecta is not effectively mixed 
into the disks, but that solid grains from the supernova likely are, thereby explaining the 
presence of SLRs hke ^°Fe in the early Solar System. 



2. Perseus 

We have written a 2-D (cylindrical) hydrodynamics code we call Perseus. Perseus (son 
of Zeus) is based heavily on the Zeus algorithms (Stone & Norman 1992). The code evolves 
the system while obeying the equations of conservation of mass, momentum and energy: 

^ + pV-v = Q (1) 

Dv , , 

p— = -Vp - pV$ (2) 

where p is the mass density, v is the velocity, p is the pressure, e is the internal energy 
density and $ is the gravitational potential (externally imposed). The Lagrangean, or 
comoving derivative D/Dt is defined as 

The pressure and energy arc related by the simple equation of state appropriate for the 
ideal gas law, p = e(7 — 1), where 7 is the adiabatic index. The term pV • v represents 
mechanical work. 
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Currently, the only gravitational potential $ used is a simple point source, representing 
a star at the center of a disk. This point mass is constrained to remain at the origin. 
Technically this violates conservation of momentum by a minute amount by excluding 
the gravitational force of the disk on the central star. As discussed in §4, the star should 
acquire a velocity ~ lO^cms"^ at the end of our simulations. In future simulations we will 
include this effect, but for the problem explored here this is completely negligible. 

The variables evolved by Perseus are set on a cyhndrical grid. The program is separated 
in two steps: the source and the transport step. The source step calculates the changes 
in velocity and energy due to sources and sinks. Using finite difference approximations, it 
evolves v and e according to 

dv 

= -Vp-pV^-V-Q (5) 



and 



Bp 

p-^^-pV-v-Q:Vv, (6) 



where Q is the tensor artificial viscosity. Detailed expressions for the artificial viscosity can 
be found in Stone & Norman (1992). 

The transport step evolves the variables according to the velocities present on the 
grid. For a given variable A, the conservation equation is solved, using finite difference 
approximat ions : 

di 



AdV = - (t Av-dS. (7) 
V Js 



The variables A advected in this way are density p, linear and angular momentum pv and 
RpVtj), and energy density e. As in the Zeus code, A on each surface element is found with 
an upwind interpolation scheme; we use second-order van Leer interpolation. 

Perseus is an explicit code and must satisfy the Courant-Friedrichs-Lewis (CFL) 
stability criterion. The amount of time advanced per timestep, essentially, must not exceed 
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the time it could take for information to cross a grid zone in the physical system. In every 
grid zone, the thermal time step Stc^ = Ax/{cs) is computed, where Ax is the size of the 
zone (smallest of the r and z dimension) and is the sound speed. Also computed are 
Str — 5r/{\vr:\) and Stz — Sz/{\vz\), where Ar and Az are the sizes of the zone in the r and 
z directions respectively. Because of artificial viscosity, a viscous time step must also be 
added for stability. For a given grid zone, the viscous time step 5tvisc = max(|(/ V ■ v/6r^)\, 
|(/ V • v/Sz"^)]) is computed, where I is a length chosen to be a 3 zone widths. The final At 
is taken to be 

At = Co {8t-^ + 8t;' + 8t-' + 8t-ZV\ (8) 

where Co is the Courant number, a safety factor, taken to be Co=0.5. To insure stability. 
At is computed over all zones, and the smallest value is kept for the next step of the 
simulation. 

Boundary conditions were implemented using ghost zones as in the Zeus code. To 
allow for supernova cjecta to flow past the disk, inflow boundary conditions were used at 
the upper boundary [z — -Zmax), and outflow boundary conditions were used at the lower 
boundary {z — ^min) and outer boundary (r = rmax)- Reflecting boundary conditions, 
were used on the inner boundary (r = rmin ^ 0) to best model the symmetry about the 
protoplanetary disk's axis. The density and velocity of gas flowing into the upper boundary 
were varied with time to match the ejecta properties (see §3). 

A more detailed description of the algorithms used in Perseus can be found in Stone & 
Norman (1992). 
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2.1. 



Additions to Zeus 



To consider the particular problem of high-velocity ejecta hitting a protoplanetary 
disk, we wrote Perseus with the following additions to the Zeus code. One minor change 
is the use of a non-uniform grid. In all of our simulations we used an orthogonal grid 
with uniform spacing in r but non-uniform spacing in the z direction. For example, in the 
canonical simulation (§3), the computational domain extends from r = 4 to 80 AU, with 
spacing Ar = 1 AU, for a total of 76 zones in r. The computational domain extends from 
z — — 50AU to -1-90 AU, but zone spacings vary with z^ from = 0.2 AU at ^ = 0, to 
t^z ~ 3 AU at the upper boundary. Grid spacings increased geometrically by 5% per zone, 
for a total of 120 zones in z. 

Another addition was the use of a radiative cooling term. The simulations bear out 
the expectation that almost all of the shocked supernova ejecta flow past the disk before 
they have time to cool significantly. Cooling is significant only where the ejecta collide with 
the dense gas of the disk itself, but there the cooling is sensitive to many unconstrained 
physical properties to do with the the chemical state of the gas, properties of dust, etc. 
To capture the gross effects of cooling (especially compression of gas near the dense disk 
gas) in a computationally simple way, we have adopted the following additional term in the 
energy equation, implemented in the source step: 



where rZe and rip are the number of protons and electrons in the gas, and A is the cooling 
function. The densities and rip are obtained simply by assuming the hydrogen gas is 
fully ionized, so rie — Up — p/lAmu- For gas temperatures above 10^ K, we take A of 
a solar-mctallicity gas from Sutherland and Dopita (1993); A typically ranges between 
10-24 erg cm^s"^ (at T = lO^K) and A = lO'^i ergcm^ s'^ (at T = 10^ K). Below lO^K 
we adopted a flat cooling function of A = lO^^^ergcm^s"^. At very low temperatures it 
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is necessary to include heating processes as well as cooling, or else the gas rapidly cools 
to unreasonable temperatures. Rather than handle transfer of radiation from the central 
star, we defined a minimum temperature below which the gas is not allowed to cool: 
Tinin = 300(r/lAU)-3/4K. Perseus uses a simple first-order, finite-difference equation to 
handle cooling. Although this method is not as precise as a predictor-corrector method, in 
§2.4 we show that it is sufficiently accurate for our purposes. 

Because Perseus is an explicit code, the implementation of a cooling term demands the 
introduction of a cooling time step to insure that the gas doesn't cool too rapidly during 
one time step, resulting in negative temperatures or other instabilities. For a radiating gas, 
the cooling timescale can be approximated by icooi ~ /csT'/nA, where ks is the Boltzmann 
constant, T is the temperature of the gas, n is the number density and A is the appropriate 
cooling function. This cooling timescale is calculated on all the grid zones where the 
temperature exceeds 10^ K, and the coohng time step 5tcooi is defined to be 0.025 times 
the shortest coohng timescale on the grid. If the smallest coohng time step is shorter that 
the previously calculated At as defined by eq. [8], then it becomes the new time step. We 
ignore zones where the temperature is below 10^ K because heating and cooling are not fully 
calculated anyway, and because these zones are always associated with very high densities 
and cool extremely rapidly, on timescales as short as hours, too rapidly to reasonably follow 
anyway. 

Finally, to follow the evolution of the ejecta gas with respect to the disk gas, a tracer 
"color density" was added. By defining a different density, the color density Pc, it is possible 
to follow the mixing of a two specific parts of a system, in this case the ejecta and the disk. 
By comparing pc to p, it is possible to know how much of the ejecta is present in a given 
zone relative to the original material. It is important to note that pc is a tracer and does 
not affect the simulation in any way. 



2.2. Sod Shock-Tube 



We have benchmarked the Perseus code against a weU-known analytic solution, the 
Sod shock tube (Sod 1978). Tests were performed to verify the validity of Perseus's results. 
It is a 1-D test, and hence was only done in the z direction, as curvature effects would 
render this test invalid in the r direction. Therefore, the gas was initially set spatially 
uniform in r. 120 zones were used in the z direction. The other initial conditions of the Sod 
shock-tube are as follows: the simulation domain is split in half and filled with a 7=1.4 gas; 
in one half {z < 0.5 cm), the gas has a pressure of 1.0 dyne cm~^ and a density of 1.0 g cm~^, 
while in the other half {z > 0.5 cm) the gas has a pressure of 0.1 dyne cm^^ and a density 
of 0.125 gcm~^. The results of the simulation and the analytical solution at t = 0.245 s 
are shown in Figure [H The slight discrepancies between the analytic and numerical results 
are attributable to numerical diffusion associated with the upwind interpolation (see Stone 
& Norman 1992), match the results of Stone & Norman (1992) almost exactly, and are 
entirely acceptable. 



2.3. Gravitational Collapse 

As a test problem involving curvature terms, we also simulated the pressure-free 
gravitational collapse of a spherical clump of gas. A uniform density gas (p = 10~'^'^gcm~^) 
was imposed everywhere within 30 AU of the star. As stated above, the only source of 
gravitational acceleration in our simulations is the central protostar, with mass M = IMq. 
The grid on which this simulation takes place has 120 zones in the z direction and 80 in the 
r direction The free-fall timescale under the gravitational potential of a 1 M0 star is 29.0 
yrs. The results of the simulation can be seen in Figure [2l After 28 years, the 30 AU clump 
has contracted to the edge of the computational volume. Spherical symmetry is maintained 
throughout as the gas is advected despite the presence of the inner boundary condition. 
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2.4. Cooling 

To test the accuracy of the coohng algorithm, a simple 2D grid of 64 zones by 64 zones 
was set up. The simulation starts with gas at T = 10^° K. The temperature of the gas is 
followed until it reaches T = 10^ K. Simulations were run varying the cooling time step 
^tcooi- As the cooling subroutine does not use a predictor-corrector method, decreasing the 
time step increases the precision. A range of cooling time steps, varying from 10 times what 
is used in the code to 0.1 times what is used in the code, were tested. Since in the range of 
T = 10^ K — lO^'' K, the cooling rate varies with temperature (according to Sutherland & 
Dopita 1993), the size of the time step should affect the time evolution of the temperature. 
This evolution is depicted in Figure [3l from which one can see that ^tcooi used in the code 
is sufficient, as using smaller time steps gives the same result. In addition, we can see that 
even the lesser precision runs give comparably good results, as the thermal time step of the 
CFL condition prevents a catastrophically rapid cooling. The precision of the cooling is 
limited by the accuracy of the cooling factors used, not the algorithm. 

2.5. "Relaxed Disk" 

Finally, we have modeled the long-term evolution of an isolated protoplanetary disk. 
To begin, a minimum-mass solar nebula disk (Hayashi et al. 1985) in Keplerian rotation 
is truncated at 30 AU. The code then runs for 2000 years, allowing the disk to find its 
equilibrium configuration under gravity from the central star (IM©), pressure and angular 
momentum. We call this the "relaxed disk" , and use it as the initial state for the runs that 
follow. To check the long term stability of the system, we allow the relaxed disk to evolve 
an extra 2000 years. This test verifies the stability of the simulated disk against numerical 
effects. In addition, using a color density, we can assess how much numerical diffusion 
occurs in the code. 
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After the extra 2000 years, the disk maintains its shape, and is deformed only at 
its lowest isodensity contour, because of the gravitational infall of the surrounding gas 
(Figure H]). Comparing this deformation to the results from the canonical run (§3), this 
is a negligible effect. Some of the surrounding gas has accreted on the disk due to the 
gravitational potential of the central star. The color density allows us to follow the location 
of the accreted gas. After 2000 years, roughly 20% of the accreted mass has found its way 
to the midplane of the disk due to the effects of numerical diffusion. Hence some numerical 
diffusion exists and must be considered in what follows. 

3. Canonical Case 

In this section, we adopt a particular set of parameters pertinent to the disk and the 
supernova, and follow the evolution of the disk and ejecta in some detail. The simulation 
begins with our relaxed disk (§2.5), seen in Figure O Its mass is about 0.00838 Mq, and 
it extends from 4AU to 40 AU, the inner parts of the disk being removed to improve 
code performance. The gas density around the disk is taken to be a uniform lOcm"'^, 
which is a typical density for an H ll region. This disk has similar characteristics to those 
found in the Orion nebula, which have been photoevaporated down to tens of AU by the 
radiation of nearby massive O stars (Johnstone, Hollenbach & Bally 1998). In setting up 
our disk, we have ignored the effects of the UV flash that accompanies the supernova, 
in which approximately 3 x 10^^ erg of high-energy ultraviolet photons are emitted over 
several days (Hamuy et al. 1988). The typical UV opacities of protoplanetary disk dust 
are k ~ lO^cm^g^^ (D'Allesio et al. 2006), so this UV energy does not penetrate below 
a column density ~ k,~^ ~ 10~^gcm~^. The gas density at the base of this layer is 
typically p ~ 10~^^gcm~^; if the gas reaches temperatures < 10^ K, tcooi will not exceed 
a few hours (§2.1). The upper layer of the disk absorbing the UV is not heated above a 
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temperature T ~ (£'uv/47rd^)mHK/A;B ~ 10^ K. Because the gas in the disk absorbs and 
then reradiates the energy it absorbs from the UV flash, we have ignored it. We have also 
neglected low-density gas structures that are likely to have surrounded the disk, including 
photoevaporative flows and bow shocks from stellar winds, as these are beyond the scope 
of this paper. It is likely that the UV flash would greatly heat this low-density gas and 
cause it to rapidly escape the disk anyway. Our "relaxed disk" initial state is a reasonable, 
simplified model of the disks seen in H ll regions before they are struck by supernova shocks. 

After a stable disk is obtained, supernova ejecta are added to the system. The canonical 
simulation assumes Mej = 20 of material was ejected isotropically by a supernova 
d — 0.3 pc away, with an explosion kinetic energy E^^ — 10^^ erg, (1 f.o.e.). This is typical of 
the mass ejected by a 25 progenitor star, as considered by Woosley & Weaver (1995), 
and although more recent models show that progenitor winds are likely to reduce the ejecta 
mass to < 10 Mq (Woosley, Heger & Weaver 2002), we retain the larger ejecta mass as a 
worst-case scenario for disk survivability. The ejecta are assumed to explode isotropically, 
but with density and velocity decreasing with time. The time dependence is taken from 
the scaling solutions of Matzner & McKee (1999); in analogy to their eq. [1], we define the 
following quantities: 



'2i^cj 
M~ 



t. = f (10) 



3Mej 



where i?* is the radius of the exploding star, taken to be 5Oi?0. The travel time from the 
supernova to the disk is computed as ttrav = c?/^*, and is typically ~ 100 years. Finally, 
expressions for the time dependence of velocity, density and pressure of the ejecta, are 
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obtained for any given time t after the shock strikes the disk: 



Vei (t) = 




4 



3 



(11) 



We acknowledge that supernova ejecta are not distributed homogeneously within the 
progenitor (Matzner & McKee 1999), nor are they ejected isotropically (Woosley, Heger 
& Weaver 2002), but more detailed modeling lies beyond the scope of this paper. 
Our assumption of homologous expansion is in any case a worst-case scenario for disk 
survivability in that the ejecta are front-loaded in a way that overestimates the ram pressure 
(C. Matzner, private communication). As our parameter study (§5) shows, density and 
velocity variations have little influence on the results. 

The incoming ejecta and the shock they create while propagating through the 
low-density gas of the H II region can be seen in Figure El When the shock reaches the disk, 
the lower-density outer edges are swept away, as the ram pressure of the ejecta is much 
higher than the gas pressure in those areas. However, the shock stalls at the higher density 
areas of the disk, as the gas pressure is higher there. A snapshot of the stalling shock can 
be seen in Figure [3 As the ejecta hit the disk, they shock and thermalize, heating the 
gas on the upper layers of the disk. This increases the pressure in that area, causing a 
reverse shock to propagate into the incoming ejecta. The reverse shock will eventually stall, 
forming a bow shock around the disk (Figures [8] and [9]). Roughly 4 months have passed 
between the initial contact and the formation of the bow shock. 

Some stripping of the low density gas at the disk's edge (> 30 AU) may occur as the 
supernova ejecta is deflected around it, due primarily to the ram pressure of the ejecta. As 
the stripped gas is removed from the top and the sides of the disk, it either is snowplowed 
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away from the disk if enough momentum has been imparted to it, or it is pushed behind 
the disk, where it can fall back onto it (Figure [TU]) . In addition to stripping the outer 
layers of the disk, the pressure of the thermalized shocked gas will compress the disk to a 
smaller size; although they do not destroy the disk, the ejecta do temporarily deform the 
disk considerably. Figure [11] shows the effect of the pressure on the disk, which has been 
reduced in thickness and has shrunk to a radius of 30 AU. The extra external pressure 
effectively aids gravity and allows the gas to orbit at a smaller radius with the same angular 
momentum. As the ejecta is deflected across the top edge of the disk, some mixing between 
the disk gas and the ejecta may occur through Kelvin-Helmholtz instabilities. Figure [T2l 
shows a close up of the disk where a Kelvin-Helmholtz roll is occurring at the boundary 
between the disk and the flowing ejecta. In addition, some ejecta mixed in with the stripped 
material under the disk might also accrete onto the disk. As time goes by and slower ejecta 
hit the disk, the ram pressure affecting the disk diminishes, and the disk slowly returns to 
its original state, recovering almost completely after 2000 years (Figure [T3l) . 

The exchange of material between the disk and the ejecta is mediated through the 
ejecta-disk interface, which in our simulations is only moderately well resolved. As discussed 
in §4, the numerical resolution will affect how well we quantify both the destruction of the 
disk and the mixing of ejecta into the disk. In the canonical run, at least, disk destruction 
and gas mixing are minimal. Although some stripping has occurred while the disk was being 
hit by the ejecta, it has lost less than 0.1% of its mass. The final disk mass, computed from 
the zones where the density is greater than 100 cm~^, remains roughly at 0.00838 M0. Some 
of the ejecta have also been mixed into the disk, but only with very low efficiency. A 30 AU 
disk sitting 0.3 pc from the supernova intercepts roughly one part in 1.7 x 10^ of the total 
ejecta from the supernova, assuming isotropic ejecta distribution. For 20 Mq of ejecta, this 
corresponds to roughly 1.18 x 10"^ M© intercepted. At the end of the simulation, we find 
only 1.48 x 10^*^ Mq of supernova ejecta was injected in the disk, for an injection efficiency 
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of about 1.3%. Some of the injected material could be attributed to numerical diffusion 
between the outer parts of the disk and the inner layers: as seen in §2.5, Perseus is diffusive 
over long periods of time. However, the distribution of the colored mass is qualitatively 
different from that obtained from a simple numerical diffusion process. Figure [T3] compares 
the percentage of colored mass within a given isodensity contour for the canonical case and 
the relaxed disk simulation of §2.5, at a time 500 years after the beginning of each of these 
simulations. From this graph, it is clear that the process that injects the supernova ejecta is 
not simply numerical diffusion, as it is much more efficient at injecting material deep within 
the disk. The post-shock pressure of the ejecta gas, 100 years after initial contact, when 
its forward progession in the disk has stalled is ~ 2pcj'i^ej/(7 + 1) = 2.8 x 10"^ dyne cm~^. 
(After 100 years, pej = 2.2 x lO^^^gcm^^ and fej = 1300 kms^^.) The shock stalls where 
the post-shock pressure is comparable to the disk pressure ~ p/c^T/m. Hence at 20 AU, 
where the temperature of the disk is T ^ 30 K, the shock stalls at the isodensity contour 
~ 1.5 X 10~^'^gcm~^. As about half of the color mass is mixed to just this depth, this is 
further evidence that the color field in the disk represents a real physical mixing. 



4. Numerical Resolution 

The results of canonical run show many similarities to related problems that have 
been studied extensively in the literature. The interaction of a supernova shock with a 
protoplanetary disk resembles the interaction of a shock with a molecular cloud, as modeled 
by Nittmann et al. (1982), Bedogni & Woodward (1990), Klein, McKee & Colella (1994; 
hereafter KMC), Mac Low et al. (1994), Xu & Stone (1995), Orlando et al. (2005) and 
Nakamura et al. (2006). Especially in Nakamura et al. (2006), the numerical resolutions 
achieved in these simulations are state-of-the-art, reaching several xlO^ zones per axis. 
In those simulations, as in our canonical run, the evolution is dominated by two physical 
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effects: the transfer of momentum to the cloud or disk; and the onset of Kelvin- Helmholtz 
(KH) instabilities that fragment and strip gas from the cloud or disk. KH instabilities 
are the most difficult aspect of either simulation to model, because there is no practical 
lower limit to the lengthscales on which KH instabilities operate (they are only suppressed 
at scales smaller than the sheared surface). Increasing the numerical resolution generally 
reveals increasingly small-scale structure at the interface between the shock and the cloud 
or disk (see Figure 1 of Mac Low et al. 1994). The numerical resolution in our canonical 
run is about 100 zones per axis; more specifically, there are about 26 zones in one disk 
radius (of 30 AU), and about 20 zones across two scale heights of the disk (one scale-height 
being about 2 AU at 20 AU). Our highest- resolution run used about 50 zones along the 
radius of the disk, and placed about 30 zones across the disk vertically. In the notation of 
KMC, then, our simulations employ about 20-30 zones per cloud radius, a factor of 3 lower 
than the resolutions of 100 zones per cloud radius argued by Nakamura et al. (2006) to be 
necessary to resolve the hydrodynamics of a shock hitting a molecular cloud. 

Higher numerical resolutions are difficult to achieve; unlike the case of a supernova 
shock with speed ~ 2000 km s~^ striking a molecular cloud with radius of 1 pc, our 
simulations deal with a shock with the same speed striking an object whose intrinsic 
lengthscale is ~ 0.1 AU. Satisfying our CFL condition requires us to use timesteps that are 
only ~ 10^ s, four orders of magnitude smaller than the timesteps needed for the case of a 
molecular cloud. This and other factors conspire to make simulations of a shock striking 
a protoplanetary disk about 100 times more computationally intensive than the case of a 
shock striking a molecular cloud. Due to the numerous lengthscales in the problem imposed 
by the star's gravity and the rotation of the disk, it is not possible to run the simulations 
at low Mach numbers and then scale the results to higher Mach numbers. We intend to 
create a parallelized version of Perseus to run on a computer cluster in the near future, but 
until then, our numerical resolution cannot match that of simulations of shocks interacting 
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with molecular clouds. This begs the question, if our resolution is not as good as has been 
achieved by others, is it good enough? 

To quantify what numerical resolutions are sufficient, we examine the physics of a 
shock interacting with a molecular cloud, and review the convergence studies of the same 
undertaken by previous authors. In the most well-known simulations (Nittmann et al. 
1982; KMC; Mac Low et al. 1994; Nakamura et al. 2006), it is assumed that a low-density 
molecular cloud with no gravity or magnetic fields is exposed to a steady shock. The 
shock collides with the cloud, producing a reverse shock that develops into a bow shock; 
a shock propagates through the cloud, passing through it in a "cloud-crushing" time tec- 
The cloud is accelerated, but as long as a velocity difference between the high-velocity gas 
and the cloud exists, KH instabilities grow that create fragments with significant velocity 
dispersions, ~ 10% of the shock speed (Nakamura et al. 2006). Cloud destruction takes 
place before the cloud is fully accelerated, and the cloud is effectively fragmented in a few 
X tec before the velocity difference diminishes. These fragments are not gravitationally 
bound to the cloud and easily escape. As long as the shock remains steady for a few x tec, 
it is inevitable that the cloud is destroyed. 

As KH instabilities are what fragment the cloud and accelerate the fragments, it is 
important to model them carefully, with numerical resolution as high as can be achieved. 
KMC stated in their abstract and throughout their paper that 100 zones per cloud 
radius were required for "accurate results" ; however, all definitions of what was meant 
by "accurate", or what were the physically relevant "results" were deferred to a future 
"Paper II". A companion paper by Mac Low et al. (1994) referred to the same Paper II 
and repeated the claim that 100 zones per axis were required. Nakamura et al. (2006), 
published this year, appears to be the Paper II that reports the relevant convergence 
study and quantifies what is meant by accurate results. Global quantities, including the 
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morphology of the cloud, its forward mass and momentum, and the velocity dispersions 
of cloud fragments, were defined and calculated at various levels of numerical resolution. 
These were then compared to the same quantities calculated using the highest achievable 
resolutions, about 500 zones per cloud radius (over 1000 zones per axis). The quantities 
slowest to converge with higher numerical resolution were the velocity dispersions, probably, 
they claim, because these quantities are so sensitive to the hydrodynamics at shocks and 
contact discontinuities where the code becomes first-order accurate only. The velocity 
dispersions converged to within 10% of the highest-resolution values only when at least 
100 zones per cloud radius were used. For this single arbitrary reason, Nakamura et al. 
(2006) claimed numerical resolutions of 100 zones per cloud radius were necessary. We note, 
however, that the other quantities to do with cloud morphology and momentum were found 
to converge much more readily; according to Figure 1 of Nakamura et al. (2006), numerical 
resolutions of only 30 zones per cloud radius are sufficient to yield values within 10% of the 
values found in the highest-resolution simulations. And although the velocity dispersions 
are not so well converged at 30 zones per cloud radius, even then the errors do not exceed 
a factor of 2. Assuming that the problem we have investigated is similar enough to that 
investigated by Nakamura et al. (2006) so that their convergence study could be applied 
to our problem, we would conclude that even our canonical run is sufficiently resolving 
relevant physical quantities, the one possible exception being the velocities of fragments 
generated by KH instabilities, where the errors could be a factor of 2. 

Of course, the problem we have investigated, a supernova shock striking a 
protoplanetary disk, is different in four very important ways from the cases considered by 
KMC, Mac Low et al. (1994) and Nakamura et al. (2006). The most important fundamental 
difference is that the disk is gravitationally bound to the central protostar. Thus, even 
if gas is accelerated to supersonic speeds ~ lOkms"''^, it is not guaranteed to escape the 
star. Second, the densities of gas in the disk, Pdisk, are significantly higher than the density 
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in the gas colliding with the disk, pej- In the notation of KMC, x — Pdisk/Pej- Because 
the disk density is not uniform, no single value of x applies, but if x is understood to 
refer to different parcels of disk gas, x would vary from lO'' to over 10^. This affects the 
magnitudes of certain variables (see, e.g.. Figure 17 of KMC regarding mix fractions), but 
also qualitatively alters the problem: the densities and pressures in the disk are so high that 
the supernova shock cannot cross through the disk, instead stalling at several scale heights 
above the disk. Unlike the case of a shock shredding a molecular cloud, the cloud-crushing 
timescale ^cc is not even a relevant quantity for our calculations. The third difference is 
that shocks cannot remain non-radiative when gas is as dense as it is near the disk. Using 
p = 10~^''gcm~^ and A = lO^^'^ergcm^s"^, tcooi is only a few hours, and shocks in the disk 
are effectively isothermal. Shocks propagating into the disk therefore stall at somewhat 
higher locations above the disk than they would have if they were adiabatic. Finally, the 
fourth fundamental difference between our simulations and those investigated in KMC, 
Mac Low et al. (1994) and Nakamura et al. (2006) is that we do not assume steady shocks. 
For supernova shocks striking protoplanetary disks about 0.3 pc away, the most intense 
effects are felt only for a time ~ 10^ years, and after only 2000 years the shock has for all 
purposes passed. There are limits, therefore, to the energy and momentum that can be 
delivered to the disk. Very much unlike the case of a steady, non-radiative shock striking 
a low-density, gravitationally unbound molecular cloud, where ultimately destruction of 
the cloud is inevitable, many factors contribute to the survivability of protoplanetary disks 
struck by supernova shocks. 

This conclusion is borne out by a resolution study wc have conducted that shows 
that the vertical momentum delivered to the disk is certainly too small to destroy it, and 
that we are not significantly underresolving the KH instabilities at the top of the disk. 
Using the parameters of our canonical case, we have conducted 6 simulations with different 
numerical resolutions. The resolutions range from truly awful, with only 8 zones in the 
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radial direction (Ar = 10 AU) and 18 zones in the vertical direction (with Az = 1 AU at 
the midplane, barely sufficient to resolve a scale height), to our canonical run (76 x 120), 
to one high-resolution run with 152 radial zones (Ar = 0.5 AU) and 240 vertical zones 
(A^ = 0.13 AU at the midplane). On an Apple G5 desktop with two 2.0-GHz processors, 
these simulations took from less than a day to 80 days to run. To test for convergence, 
we calculated several global quanities Q, including: the density-weighted cloud radius, a; 
the density- weighted cloud thickness, c; the density- weighted vertical velocity, (f^); the 
density- weighted velocity dispersion in r, 5vr] the density- weighted velocity dispersion in z, 
{vz)', as well as the mass of ejecta injected into the disk, Mjnj. Except for the last quantity, 
these are defined exactly as in Nakamura et al. (2006), but using a density threshold 
corresponding to 100 cm~^. Each global quantity was measured at a time 500 years into 
each simulation. We define each global quantity Q as a function of numerical resolution n, 
where n is the geometric mean of the number of zones along each axis, which ranges from 
12 to 191. To compare to the resolutions of KMC, one must divide this number by about 3 
to get the number of zones per "cloud radius" (two scale heights at 20 AU) in the vertical 
direction, and divide by about 2 to get the number of zones per cloud radius in the radial 
direction. The convergence is measured by computing \Q{n) — (5(nmax)| /Q('^max), where 
^max = 191 corresponds to our highest resolution case. In Figure [15] we plot each quantity 
Q{n) as a function of resolution n (except (f^)). All of the quantities have converged 
to within 10%, the criterion imposed by Nakamura et al. (2006) as signifying adequate 
convergence. It is significant that 5vr has converged to within 10%, because this is the 
quantity relevant to disk destruction by KH instabilities. Material is stripped from the 
disk only if supersonic gas streaming radially above the top of the disk can generate KH 
instabilities and fragments of gas that can then be accelerated radially to escape velocities. 
If we were underresolving this layer significantly, one would expect large differences in 5vr 
as the resolution was increased, but instead this quantity has converged. Higher-resolution 
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simulations are likely to reveal smaller-scale KH instabilities and perhaps more stripping of 
the top of the disk, but not an order of mangitude more. 

The convergence of {vz) with resolution is handled differently because unlike the other 
quantities, {vz) can vanish at certain times. The disk absorbs the momentum of the ejecta 
and is pushed downward, but unlike the case of an isolated molecular cloud, the disk feels 
a restoring force from the gravity of the central star. The disk then undergoes damped 
vertical oscillations about the origin as it collides with incoming ejecta at lower and lower 
speeds. This behavior is illustrated by the time-dependence of (f^), shown in Figure [16] 
for two numerical resolutions, our canonical run {n = 95) and our highest-resolution run 
(n = 191). Figure [16] shows that the vertical velocity of the disk oscillates about zero, but 
with an amplitude ~ 0.1 km s^^. The time-average of this amplitude can be quantified by 



\^{vz) — < Vz , where the bar represents an average over time; the result is 825 cms~^ 

for the highest-resolution run and is only 2% smaller for the canonical resolution. The 
difference between the two runs is generally much smaller than this; except for a few times 
around t = 150 yr, and t = 300 yr, when the discrepancies approach 30%, the agreement 
between the two resolutions is within 10%. The time-averaged dispersion of the amplitude 
of the difference (defined as above for (vz) itself) is only 12.0 cm s~^, which is only 1.5% of 
the value for {vz) itself. Taking a time average of \{vz)g^ — {vz)igi\ / |(^z)i9il yields 8.7%. 
We therefore claim convergence at about the 10% level for (vz) as well. 

Using these velocities, we also note here that the neglect of the star's motion is entirely 
justified. The amplitude of (vz) is entirely understandable as reflecting the momentum 
delivered to the disk by the supernova ejecta, which is ~ 20 Mq {nR'^^^^^/AncP) Vej ~ 
10^'^ km s~^, and which should yield a disk velocity ~ 0.1 km s^^. The period of 
oscillation is about 150 years, which is consistent with most of this momentum being 
delivered to the outer reaches of the disk from 25 to 30 AU where the orbital periods are 
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125 to 165 years. These velocities are completely unaffected by the neglected velocity of the 
central star, whose mass is 120 times greater than the disk's mass. If the central star, with 
mass ~ 1 Mr,,, had been allowed to absorb the ejecta's momentum, it would only move at 
~ lOOcms"^ and be displaced at most 0.4 AU after 2000 years. This neglected velocity, is 
much smaller than all other relevant velocities in the problem, including \ {vz) \ ~ 800cms~^, 
as well as the escape velocities (~ lOkms"^), the velocities of gas flowing over the disk 
(~ lO^kms"^), and of course the shock speeds (~ lO^kms"^). 

Our analysis shows that wc have reached adequate convergence with our canonical 
numerical resolution {n = 95) . We observe KH instabilities in all of our simulations (except 
n — 12), and we see the role they play in stripping the disk and mixing eject a gas into it. 
We are therefore confident that we are adequately resolving these hydrodynamic features; 
nevertheless, we now consider a worst-case scenario in which we KH instabilities can strip 
the disk with 100% efficiency where they act, and ask how much mass the disk could 
possibly lose under such conditions. 

Supernova ejecta that has passed through the bow shock and strikes the disk necessarily 
stalls where the gas pressure in the disk exceeds the ram pressure of the ejecta. Below this 
level, the momentum of the ejecta is transferred not as a shock but as a pressure (sound) 
wave. Gas motions below this level are subsonic. Note that this is drastically different from 
the case of an isolated molecular cloud as studied by KMC and others; the high pressure in 
the disk is maintained only because of the gravitational pull of the central star. 

The location where the incoming ejecta stall is easily found. Assuming the vertical 
isothermal minimum-mass solar nebula disk of Hayashi et al. (1985), the gas density varies 
as p(r, z) = 1.4 x 10~^ (r/1 AU)~^^/^ exp(— ;s^/2i7^) g cm~^, where H = Cg/fi, Cg is the sound 
speed and Q is the Keplerian orbital frequency. Using the maximum density and velocity 
of the incoming ejecta (pej = 1-2 x 10~^°gcm~^ and V^j = 2200 km s"^), the ram pressure 
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of the shock striking the disk does not exceed pram = PejKj/^ — 1-5 x 10~^dynecm~^ (the 
factor of 1/4 arises because the gas must pass through the bow shock before it strikes the 
disk). At 10 AU the pressure in the disk, pc^, exceeds the ram pressure at z = 2.7 H, and at 
20 AU the ejecta stall eit z — 1.7H; the gas densities at these locations are 10~^^gcm~^. 
At later times, the ejecta stall even higher above the disk, because pram oc (cf. eq. [11]). 
For example, at t = 100 yr, the ram pressure drops below 1 x 10"^ dyne cm~^, and the ejecta 
stall above z = 3.6H (10 AU) and z = 2.9H (20 AU). 

The column density above a height z in a. vertically isothermal disk is easily found to 
be S(> z) ~ p{z)H'^/z — p{z) / {fl"^ z) . Integrating over radius, the total amount of disk gas 
that ever comes into contact with ejecta is (approximating z — 2H): 



Using a disk radius Rd — 30 AU, the maximum amount of disk gas that is actually 
exposed to a shock at any time is only 1.5 x 10~^ Mq, or 0.2% of the disk mass. This 
fraction decreases with time as Pram oc (eq. [11]); the integral over time of Pram is 
Pram(^ = 0) X itrav/4. The ram pressure drops so quickly, that effectively ejecta interact 
with this uppermost 0.2% of the disk mass only for about 30 years. This is equivalent to 
one orbital timescale at 10 AU, so the amount of disk gas that is able to mix or otherwise 
interact with the ejecta hitting the upper layers of the disk is very small, probably a few 
percent at most. As for KH instabilities, they are initiated when the Richardson number 
drops below a critical value, when 



where g — —D.'^z is the vertical gravitational acceleration, Q is the Keplerian orbital 
frequency, and (dU/Oz) is the velocity gradient at the top of the disk. Below the stall 
point, all gas motions are subsonic and the velocity gradient would have to be exceptionally 
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steep, with an unreasonably thin shear layer thickness, ^i//10, to initiate KH instabilities. 
Mixing of ejecta into the disk is quite effective above where the shock stalls, as illustrated 
by Figure dH it is in these same layers (experiencing supersonic velocities) that we expect 
that KH instabilities to occur, but again < 1% of the disk mass can be expected to interact 
with these layers. 

To summarize, our numerical simulations are run at a lower simulation (by a factor 
of about 3) than has been claimed necessary to study the interaction of steady shocks 
with gravitationally unbound molecular clouds, but the drastically different physics of the 
problem studied here as allowed us to achieve numerical convergence and allowed us to reach 
meaningful conclusions. Our global quantities have converged to within 10%, the same 
criterion used by Nakamura et al. (2006) to claim convergence. The problem is so different 
because the disk is tightly gravitationally bound to the star and the supernova shock is 
of finite duration. The high pressure in the disk makes the concept of a cloud-crushing 
time meaningless, because the ejecta stall before they drive through even 1% of the disk 
gas. Rather than a sharp interface between the ejecta and the disk, the two interact via 
sound waves within the disk, which entails smoother gradients. While we do resolve KH 
instabilities in this interface, we allow that we may be underresolving this layer; but even 
if we are, this will not affect our conclusions regarding the disk survival or the amount of 
gas mixed into the disk. This is because we already find that mass is stripped from the 
disk and ejecta are mixed into the disk very effectively (see Figure UM above the layer 
where the ejecta stall, and below this layer mixing is much less efficient and all the gas is 
subsonic and bound to the star. It is inevitable that mass loss and mixing of ejecta should 
be only at the ~ 1% level. Similar studies using higher numerical resolutions are likely to 
reveal more detailed structures at the disk-ejecta interface, but it is doubtful that more 
than a few percent of the disk mass can be mixed-in ejecta, and it is even more doubtful 
that even 1% of the disk mass can be lost. We therefore have sufficient confidence in our 
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canonical resolution to use it to test the effects of varying parameters on gas mixing and 
disk destruction. 

5. Parameter Study 
5.1. Distance 

Various parameters were changed from the canonical case to study their effect on the 
survival of the disk and the injection efficiency of ejecta, including: the distance between the 
supernova and the disk, d; the explosion energy of the supernova, E'ej; and the mass of gas 
in the disk, Mjisk- In all these scenarios, the resolution stayed the same as in the canonical 
case. The first parameter studied was the distance between the supernova and the disk. 
Prom the canonical distance of 0.3 pc, the disk was moved to 0.5 pc and 0.1 pc. The main 
effect of this change is to vary the density of the ejecta hitting the disk (see eq. [11]). If the 
disk is closer, the gaseous ejecta is less diluted as it hits the disk. Hence these simulations 
are essentially equivalent to simulating a denser or a more tenuous clump of gas hitting the 
disk in an non-homogeneous supernova explosion. The results of these simulations can be 
seen in Table 2. The "% injected" column gives the percentage of the ejecta intercepted by 
the disk [with an assumed cross-section of 7r(30 AU)^] that was actually mixed into the disk. 
The third column gives the estimated ^^Al/^^Al ratio that one would expect in the disk 
if the SLRs were dehvered in the gas phase. This quantity was calculated using a disk 
chemical composition taken from Lodders (2003), and the ejecta isotopic composition from 
a 25 Mq supernova taken from Woosley & Weaver (1995), which ejects M — 1.27 x 10~'^Mq 
of ^^Al. Although the injection efficiency increases for denser ejecta, and the geometric 
dilution decreases for a closer supernova, gas-phase injection of ejecta into a disk at 0.1 pc 
cannot explain the SLR ratios in meteorites. The ^^Al/^^Al ratio is off by roughly an order 
of magnitude from the measured value of 5 x 10~^ (e.g., MacPherson et al. 1995). Stripping 
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was more important with denser ejecta {d — 0.1 pc), although still negligible compared to 
the mass of the disk; only 0.7% of the disk mass was lost. 

5.2. Explosion Energy 

We next varied the explosion energy, which defines the velocity at which the ejecta 
travel. The explosion energy was changed from 1 f.o.e. to 0.25 and 4 f.o.e., effectively 
modifying the ejecta velocity from 2200 km/s to 1100 km/s and 4400 km/s, respectively. 
The results of the simulations can be seen in Table 3. Slower ejecta thermalizes to a lower 
temperature, and does not form such a strong reverse shock. Therefore, slower ejecta 
is injected at a slightly higher efficiency into a disk. Primarily, though, the results are 
insensitive to the velocity of the incoming supernova ejecta. 

5.3. Disk Mass 

The final parameter varied was the mass of the disk. Prom these simulations, the mass 
of the the minimum mass disk used in the canonical simulation was increased by a factor 
of 10, and decreased by a factor of 10. The results of the simulations can be seen in Table 
4. Increasing the mass by a factor of 10 shghtly increases, but this could be due to the fact 
that the disk does not get compressed as much as the canonical disk (it has a higher density 
and pressure at each radius). Hence the disk has a larger surface to intercept the ejecta 
(the calculation for injection efficiency assumes a radius of 30 AU). Reducing the mass by 
a factor of 10 increases the efficiency. As the gas density in the disk is less, the pressure 
is less, and hence the ejecta is able to get closer to the midplane, increasing the amount 
injected. 
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In this paper, we have described a 2-D cylindrical hydrodynamics code we wrote, 
Perseus, and the results from the application of this code to the problem of the interaction 
of supernova shocks with protoplanetary disks. A main conclusion of this paper is that disks 
are not destroyed by a nearby supernova, even one as close as 0.1 pc. The robustness of 
the disks is a fundamentally new result that differs from previous 1-D analytical estimates 
(Chevalier 2000) and numerical simulations (Ouellette et al. 2005). In those simulations, in 
which gas could not be deflected around the disk, the full momentum of the supernova ejecta 
was transferred directly to each annulus of gas in the disk. Chevalier (2000) had estimated 
that disk annuli would be stripped away from the disk wherever Mgj V"ej /47r(i^ > SdKsc, 
where Ed is the surface density of the disk [Ed = 1700 (r/1 AU)^'^/^ gcm~^ for a minimum 
mass disk; Hayashi et al. 1985), and V^sc is the escape velocity at the radius of the annulus. 
In the geometry considered here, the momentum is applied at right angles to the disk 
rotation, so fesc can be replaced with the Keplerian orbital velocity, as the total kinetic 
energy would then be sufficient for escape. Also, integrating the momentum transfer over 
time (eq. [11]), we find V^j = 3t'j,/4. Therefore, using the criterion of Chevalier (2000), and 
considering the parameters of the canonical case but with d = 0.1 pc, the disk should have 
been destroyed everywhere outside of 30.2 AU, representing a loss of 13% of the mass of a 
40 AU radius disk. Comparable conclusions were reached by Ouellette et al. (2005). 

In contrast, as these 2-D simulations show, the disk becomes surrounded by high- 
pressure shocked gas that cushions the disk and defiects ejecta around the disk. This 
high-pressure gas has many effects. First, the bow shock deviates the gas, making part 
of the ejecta that would have normally hit the disk flow around it. From Figure fTTl by 
following the velocity vectors, it is possible to estimate that the gas initially on trajectories 
withr > 20 AU will be deflected by > 14° after passing through the bow showk, and will 
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miss the disk. For a disk 30 AU in size, this represents a reduction in the mass flux hitting 
by 45%; more thorough calculations give a reduction of ki 50%. Second, the bow shock 
reduces the forward velocity of the gas that does hit the disk. Gas deviated sideways about 
14°, will have lost more than 10% of its forward velocity upon reaching the disk. These 
two effects combined conspire to reduce the amount of momentum hitting the disk by 55% 
overall. By virtue of the smaller escape velocity and the lower disk surface density, gas at 
the disk edges is most vulnerable to loss by the momentum of the shock, but it at the disk 
edges that the momentum of the supernova shock is most sharply reduced. Because of 
the loss of momentum, the disk in the previous paragraph could survive out to a radius of 
about 45 AU. 

A third, significant effect of the surrounding high-pressure shocked gas, though, is 
its ability to shrink the disk to a smaller radius. The pressure in the post-shock gas is 
~ 2pejfej/(7 + 1) = 4.4 X lO^^dynecm"^, so the average pressure gradient in the disk 
between about 30 and 35 AU is fa 1.9 x 10~^^ dynecm"^. This is to be compared to 
the gravitational force per volume at 35 AU, pg = 4.8 x 10~^^ dyne cm"^ (at 35 AU, 
p ~ 1.0 X 10~^^ in the canonical disk.) The pressure of the shocked gas enhances the inward 
gravitational force by a significant amount, causing gas of a given angular momentum 
to orbit at a smaller radius than it would if in pure Keplerian rotation. When this high 
pressure is relieved after the supernova shock has passed, the disk is restored to Keplerian 
rotation and expands to its original size. While the shock is strongest, the high-pressure gas 
forces a protoplanetary disk to orbit at a reduced size, ~ 30 AU, where it is invulnerable to 
being stripped by direct transfers of momentum. Because of these combined effects of the 
cushion of high-pressure shocked gas surrounding the disk — reduction in ejecta momentum 
and squeezing of the disk — protoplanetary disks even 0.1 pc from the supernova lose < 1% 
of their mass. 
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Destruction of the disk, if it occurs at all, is due to stripping of the low-density upper 
layers of the disk by Kelvin-Helmholtz (KH) instabilities. We observed KH instabilities in 
all of our simulations (except n = 12), and we observe their role in stripping gas from the 
disk and mixing supernova ejecta into the disk (e.g.. Figure [T2!) . Our canonical numerical 
resolution [n = 95) and our highest-resolution simulation {n = 191), corresponding to 
effectively 20-30 zones per cloud radius in the terminology of KMC, are just adequate to 
provide convergence at the 10% level, as described in §4. We are confident we are capturing 
the relevant physics in our simulations, but we have shown that even if KH instabilities 
are considerably more effective than we are modeling, that no more than about 1% of 
the disk mass could ever be affected by KH instabilities. This is because the supernova 
shock stalls where the ram pressure is balanced by the pressure in the disk, and for typical 
protoplanetary disk conditions, this occurs several scale heights above the midplane. It is 
unlikely that higher-resolution simulations would observe loss of more than ~ 1% of the 
disk mass. We observed that the ratio of injected mass to disk mass was typically ~ 1% as 
well, for similar reasons. We stipulate that mixing of ejecta into the disk is more subtle 
than stripping of disk mass, but given the limited ability of the supernova ejecta to enter 
the disk, we find it doubtful that higher-resolution simulations would increase by more than 
a few the amount of gas-phase radionuclides injected into the disk. Therefore, while disks 
like those observed in the Orion Nebula (McCaughrean & O'Dell 1995) should survive the 
explosions of the massive stars in their vicinity, and while these disks would then contain 
some fraction of the supernova's gas-phase ejecta, they would not retain more than a small 
fraction (~ 1%) of the gaseous ejecta actually intercepted by the disk. If SLRs like ^^Al are 
in the gas phase of the supernova (as modeled here), they will not be injected into the disk 
in quantities large enough to explain the observed meteoritic ratios, failing by 1-2 orders of 
magnitude. 

Of course, the SLRs inferred from meteorites, e.g., ^°Fe and ^^Al, would not be detected 
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if they were not refractory elements. These elements should condense out of the supernova 
ejecta as dust grains before colliding with a disk (Ebel & Grossman 2001). Colgan et al. 
(1994) observed the production of dust at the temperature at which FeS condenses, 640 
days after SN 1987A, suggesting that the Fe and other refractory elements should condense 
out of the cooling supernova ejecta in less than a few years. (The supernova ejecta is 
actually quite cool because of the adiabatic expansion of the gas.) As the travel times from 
the supernova to the disks in our simulations are typically 20 — 500 years, SLRs can be 
expected to be sequestered into dust grains condensed from the supernova before striking a 
disk. 

Dust grains will be injected into the disk much more effectively than gas-phase ejecta. 
When the ejecta gas and dust, moving together at the same speed, encounter the bow 
shock, the gas is almost instantaneously deflected around the disk, but the dust grains will 
continue forward by virtue of their momentum. The dust grains will be slowed only as fast 
as drag forces can act. The drag force F on a highly supersonic particle is F ^ vra^ pg Ai>^, 
where a is the dust radius, pg is the gas density, and Av is the velocity difference between 
the gas and the dust. Assuming the dust grains are spherical with internal density Ps, the 
resultant acceleration is dv/dt— — (3pgAi>^)/(4psa)- Immediately after passage through the 
bow shock, the gas velocity has dropped to 1/4 of the ejecta velocity, so Av (3/4) Wej- 
Integrating the acceleration, we find the time ti/2 for the dust to lose half its initial velocity: 

Measurements of SiC grains with isotopic ratios indicative of formation in supernova ejecta 
reveal typical radii of a ~ 0.5 pm (Amari et al. 1994; Hoppe et al. 2000). Assuming 
similar values for all supernova grains, and an internal density Ps — 2.5gcm~^, and using 
the maximum typical gas density in the region between the bow shock and the disk, 
Pg ~ 5 X 10~^''gcm~^ , we find a minimum dust stopping time ti/2 ~ 2 x 10^ s. In that 
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time, the dust will have travelled about 300 AU. As the bow shock hes about 20 AU from 
the disk, the dust will encounter the protoplanetary disk well before travelling this distance, 
and we conclude that the dust the size of typical supernova SiC grains is not deflected 
around the disk. We estimate that nearly all the dust in the ejecta intercepted by the 
disk will be injected into the disk. With nearly 100% injection efficiency, the abundances 
of ^^Al and ^''Fe in a disk 0.15 pc from a supernova would be ^^A1/^''A1 = 6.8 x 10~^ 
and *'°Fe/^^Fe = 4.8 x 10"'' (using the yields from Woosley & Weaver 1995). These 
values compare quite favorably to the meteoritic ratios (^^Al/^'^'Al = 5.0 x 10"^ and 
60pg/56pg ^ 3-7x 10-^; MacPherson et al. 1995, Tachibana & Huss 2003), and we conclude 
that injection of SLRs into an already formed protoplanetary disk by a nearby supernova 
is a viable mechanism for delivering radionuclides to the early Solar System, provided the 
SLRs have condensed into dust. In future work we will present numerical simulations of 
this process (Ouellette et al. 2007, in preparation). 

We thank an anonymous referee for two very thorough reviews that significantly 
improved the manuscript. We also thank Chris Matzner for helpful discussions. 
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Table 1 Mass injected 



# of zones {r x z) % injected 



8x18 


0.83 


16x30 


0.77 


30x54 


0.96 


39x74 


1.28 


60x88 


1.31 


76 xl20 


1.26 


152x240 


1.25 



Table 2 Effect of Distance 

d % injected ^6^1/27^1 

0.1 pc 4.3 6.4 X 10"^ 

0.3 pc 1.3 2.2 X 10"^ 

0.5 pc 1.0 5.6 X 10-^ 



Table 3 Effect of Explosion Energy 

E^^ % injected ^6^1/27^1 

4.0 f.o.e. 1.0 1.7 X 10"^ 

1.0 f.o.e. 1.3 2.2 X 10"^ 

0.25 f.o.e. 1.7 2.8 x 10"^ 
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Table 4 Effect of Disk Mass 

Mdisk % injected ^6^1/27^1 

0.084 Mq 1.4 2.3 X 10"^ 

0.0084 M0 1.3 2.2 X 10"^ 

0.00084 Mr. 2.2 3.6 x 10"^ 
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Fig. 1. — Sod shock-tube problem benchmark. The squares are the results of the simulation 
using Perseus, and the solid hne is the analytical solution from Hawley et al. (1984). 
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Fig. 2. — Pressure-free collapse of a 30 AU, uniform-density clump of gas, as simulated by 
Perseus. Spherical symmetry is maintained despite the cylindrical geometry and the inner 
boundary condition at r = 2 AU. 
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Fig. 3. — (a) Gas temperature evolution using various cooling time steps {tc is the time step 
normally used in the code) (b) Close-up of the time interval 31-33 years. Convergence is 
achieved using and higher resolution timesteps. 



-45 - 



a) 40 



ZD 

< 




10 20 



30 40 
R (AU) 



50 60 



b) 40 



< 




10 20 30 40 
R (AU) 



50 60 



Fig. 4. — (a) Isodensity contours of an equilibrium ("relaxed") protoplanctary disk. Con- 
tours are spaced a factor of 10 apart, with the outermost contour representing a density 
10~^°gcm~^. The rotating disk has already evolved for 2000 years and is stable; this is the 
configuration used as the initial state for our subsequent runs, (b) The "relaxed" disk, after 
an additional 2000 years of evolution. While some slight deformation of the lowest-density 
contours is seen, attributable to gravitational infall of surrounding gas, the disk is stable and 
non-evolving over the spans of time relevant to supernova shock passage, ~ 2000 years. 
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Fig. 5. — Isodensity contours of the relaxed disk, just before impact of the supernova ejecta. 
Contours are spaced a factor of 10 apart, with the outermost contour representing a density 
10~^^ gcm~^. 
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Fig. 6. — Protoplanetary disk immediately prior to impact by the supernova shock. Isoden- 
sity contours are as in Figure HI Arrows represent gas velocities. The supernova ejecta are 
travelling through the Hll region toward the disk at about 2200 km s~^. 
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Fig. 7. — Protoplanetary disk 0.05 years after being first hit by supernova ejecta. As the 
supernova shock sweeps around the disk edge, it snowplows the low-density disk material 
with it, but the shock stalls in the high-density gas in the disk proper. Isodensity contours 
and velocity vectors as in Figure [51 
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Fig. 8. — Protoplanetary disk 0.1 years after first being hit by supernova ejecta. As the 
pressure increases on the side of the disk facing the ejecta, a reverse shock forms, visible as 
the outermost (dashed) isodensity contour. Isodensity contours and velocity vectors as in 
Figure [51 
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Fig. 9. — Protoplanetary disk 0.3 years after first being hit by supernova ejecta. The reverse 
shock, visible as the outermost (dashed) contour, has stalled and formed a bow shock. The 
bow shock deflects incoming gas around the disk, which is effectively protected in a high- 
pressure "bubble" of gas. Isodensity contours and velocity vectors as in Figure [51 
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Fig. 10. — Protoplanetary disk 4 years after first being hit by supernova ejecta. Gas is being 
stripped from the disk (e.g., the clump between R = 35 and 45 AU). Gas stripped from the 
top of the disk either is entrained in the fiow of ejecta and escapes the simulation domain, 
or fiows under the disk and falls back onto it. Isodensity contours and velocity vectors as in 
Figure [51 
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Fig. 11. — Protoplanetary disk 50 years after first being hit by supernova ejecta. The disk is 
substantially deformed by the high pressures in the surrounding shocked gas. The pressures 
compress the disk in the z direction, and also effectively aid gravity in the r direction, 
allowing the gas to orbit at smaller radii with the same angular momentum. Isodensity 
contours and velocity vectors as in Figure [51 
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Fig. 12. — (a) Protoplanetary disk 400 years after first being hit by supernova ejecta. At 
this instant mass is being stripped off the top of the disk by a Kelvin-Helmholtz instabihty, 
seen in detail in (b). Isodensity contours and velocity vectors as in Figure [5l 
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Fig. 13. — (a) Protoplanetary disk prior to impact by supernova ejecta (same as the relaxed 
disk of FigureHj), and (b) 2000 years after first being struck by supernova ejecta. This "before 
and after" picture of the disk illustrates how the disk recovers almost completely from the 
shock of a nearby supernova. Isodensity contours and velocity vectors as in Figure [51 
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Fig. 14. — Ratio of color mass to total mass within a given isodensity contour (abscissa). 
The dashed line represents the mass ratio after allowing the relaxed disk to evolve for 2000 
years in the absence of a supernova; the solid line represents the mass ratio after 2000 years 
of interaction with the supernova shock (our canonical simulation). Supernova ejecta is 
injected very effectively up to densities where the shock would stall (~ 10~^^gcm~^), much 
more effectively than can be accounted for by numerical diffusion alone. 
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Fig. 15. — Convergence properties of selected global variables. The global variables are: the 
mass-weighted radius of the cloud, a; the mass-weighted cloud thickness, c; the dispersions 
in the mass- weighted radial {Svr) and vertical {{vz)) velocities; and the mass of ejecta gas 
injected into the disk, Minj. The quantities are calculated at a time t — 500 years, but using 
6 different numerical resolutions, n = 12, 22, 40, 54, 98 and 191. The deviation of each global 
quantity Q from the highest-resolution value Qigi is plotted against numerical resolution n. 
For our canonical simulation {n — 98), all quantities have converged at about the 10% level. 
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Fig. 16. — Density-weighted velocity along the z axis using the highest numerical simulation 
n — 191 (solid hne) and the canonical resolution n = 98 (dotted hne). The difference 
between them is plotted as the dashed line. After absorbing the initial impulse of downward 
momentum from the supernova ejecta, the disk oscillates vertically about the position of the 
central protostar with a period ~ 150 years, characteristic of the most affected gas at about 
30 AU. 



